#--------------------------------------------------------------
# Script for replicating Camp-Level results
# Egerod and Junk (2022)
# 
# The output of this script will correspond to Figure 3 in the article
#
#--------------------------------------------------------------


# load some packages
library(haven); library(plyr); library(dplyr)
library(sp); library(spdep);  library(splm); library(igraph) # load packages
library(spatialprobit)
library(stargazer)
library(ggplot2)

# set working directory to script location
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#load data
df <- read_dta("CampData.dta")

#extract used data
used.dat <- dplyr::select(df, relsharestratcoalcamp, relstaff_resources_camp, relsharebusinessactorscamp_all,
                          relsharepublicactorscamp_all, shareactorscamp_all, diffhhi_campactorRM,
                          prochange_survey, posamepos_actor_survey,success, average_no_statements,
                          actorid, policyid_all, country, position_survey_plus)

# create spatial weights matrix
used.dat$position_issue <- paste(used.dat$policyid_all, used.dat$position_survey_plus, used.dat$country, sep = "_")

camp_lvl <- subset(used.dat, duplicated(position_issue) == F)
camp_lvl <- na.omit(camp_lvl)


camp_lvl$actor.issue <- paste(camp_lvl$actorid, camp_lvl$policyid_all, camp_lvl$country, sep = "-")

spat.position <- t(combn(camp_lvl$position_survey_plus, m = 2))
spat.issue <- t(combn(camp_lvl$policyid_all, m = 2))
spat.actor <- t(combn(camp_lvl$actor.issue, m = 2))
spat.cntry <- t(combn(camp_lvl$country, m = 2))
spat <- data.frame(spat.position, spat.issue, spat.actor, spat.cntry)

names(spat) <- c("position_i", "position_j", "issue_i", "issue_j", "actor_i", "actor_j", "country_i", "country_j")
spat$actor_i <- as.character(spat$actor_i)
spat$actor_j <- as.character(spat$actor_j)


connect <- ifelse(spat$position_i == spat$position_j & spat$actor_i != spat$actor_j
                  & spat$issue_i == spat$issue_j & spat$country_i == spat$country_j, 1, 0)

connect <- as.matrix(connect)

#for camp level estimates
spat.position <- t(combn(camp_lvl$position_survey_plus, m = 2))
spat.issue <- t(combn(camp_lvl$policyid_all, m = 2))
spat.cntry <- t(combn(camp_lvl$country, m = 2))

spat <- data.frame(spat.position, spat.issue, spat.cntry)

names(spat) <- c("position_i", "position_j", "issue_i", "issue_j", "country_i", "country_j")

connect <- ifelse(spat$position_i != spat$position_j
                  & spat$issue_i == spat$issue_j & spat$country_i == spat$country_j, 1, 0)

connect <- as.matrix(connect)

spat.CampId<- t(combn(camp_lvl$position_issue, m = 2))

w.edge <- data.frame(spat.CampId, connect)
w.edge[,1] <- as.character(w.edge[,1])
w.edge[,2]<- as.character(w.edge[,2])

w.mat <- as.matrix(w.edge)
w.mat[,1] <- as.character(w.mat[,1])
w.mat[,2]<- as.character(w.mat[,2])

w.mat <- as.matrix(w.mat)

w.mat2 <- graph.edgelist(w.mat[, 1:2], directed = F)
E(w.mat2)$weight <- as.numeric(as.character(w.edge$connect))

panel.weight <- get.adjacency(w.mat2, attr='weight')

slx <- as.matrix(panel.weight)
panel.weight <- mat2listw(slx)

weights_matrix <- Matrix(slx, sparse = TRUE)

panel.weight <- as(as_dgRMatrix_listw(panel.weight), "CsparseMatrix")

camp_lvl <- as.data.frame(camp_lvl)

# reverse diversity variable
camp_lvl$rev_diffhhi <- camp_lvl$diffhhi_campactorRM*-1


# estimate
set.seed(182)

mod1 <-sarprobit(success ~ relsharestratcoalcamp + 
                   relstaff_resources_camp + 
                   rev_diffhhi +
                   shareactorscamp_all +
                   posamepos_actor_survey + 
                   prochange_survey + 
                   average_no_statements, #+
                 #factor(country),
                 data = camp_lvl, W = panel.weight, thinning = 20,
                 showProgress = T, burn.in = 3500)#max burnin w/ country FE = 2000

sum_camp <- summary(mod1)
camp_se <- sum_camp[,2]
camp_coef <- sum_camp[,1]

rho_sd <- sum_camp[9,2]

prob_mod<-glm(success ~ relsharestratcoalcamp + 
                relstaff_resources_camp + 
                rev_diffhhi +
                shareactorscamp_all +
                posamepos_actor_survey + prochange_survey + 
                average_no_statements ,#+
              #factor(country),
              data = camp_lvl, 
              family = binomial(link = "probit"))

#-----------
# Table for the appendix

stargazer(prob_mod, coef= list(camp_coef),
          se = list(camp_se),
          covariate.labels = c("Strategy Coalitions",
                               "# Staff",
                               "HHI",
                               "Camp Size",
                               "Public Support",
                               "Status Quo",
                               "Salience"),
          title = "Bayesian SAR Probit Coefficients:  Opposing Camps",
          dep.var.labels = "Preference Attainment",
          add.lines = list(c("Spatial Parameter", round(mod1$rho, digits = 2)),
                           c(" ", round(rho_sd, digits = 2))),
          out = "camp_bayes_regressions.html",
          type = "html")


#-------------------
# calculate direct + indirect effects for graph

indirect <- as.data.frame(mod1$indirect)
apply(indirect, 2, quantile, probs = c(.025, 0.05, 0.5, 0.95, 0.975))

indirect$rev_diffhhi <- indirect$rev_diffhhi*sd(camp_lvl$rev_diffhhi)
indirect$shareactorscamp_all <- indirect$shareactorscamp_all*sd(camp_lvl$shareactorscamp_all)


direct <- as.data.frame(mod1$direct)
apply(direct, 2, quantile, probs = c(.025, 0.05, 0.5, 0.95, 0.975))

direct$rev_diffhhi <- direct$rev_diffhhi*sd(camp_lvl$rev_diffhhi)
direct$shareactorscamp_all <- direct$shareactorscamp_all*sd(camp_lvl$shareactorscamp_all)


total <- mod1$total
apply(total, 2, quantile, probs = c(.025, 0.05, 0.5, 0.95, 0.975))

spat_fx <- rbind(cbind(apply(direct, 2, quantile, probs = c(0.5)),
                       t(apply(direct, 2, quantile, probs = c(0.025, 0.975))),
                       t(apply(direct, 2, quantile, probs = c(0.05, 0.95)))),
                 cbind(apply(indirect, 2, quantile, probs = c(0.5)),
                       t(apply(indirect, 2, quantile, probs = c(0.025, 0.975))),
                       t(apply(indirect, 2, quantile, probs = c(0.05, 0.95)))))
spat_fx <- as.data.frame(spat_fx) 
spat_fx$effect <- c(rep("A: Direct Effects", 7),
                    rep("B: Indirect Effects", 7))
spat_fx$var <- rownames(spat_fx)

names(spat_fx)[1:5] <- c("PE", "lwr95", "upr95",
                         "lwr90", "upr90") 

spat_fx$var_num <- 1:7

ann_text_new <- data.frame(PE = c(0.3,0), var = spat_fx$var[c(3,4)],
                       var_num = c(4.15,4),ordering = c(2.09,2.09), 
                       lab = c(paste("rho==", round(mod1$rho,2),sep = ""), NA), 
                       effect = c("A: Direct Effects", 
                                  "B: Indirect Effects"))


# create graph with only the variables of interest

small_plot <- spat_fx %>% filter(var == "rev_diffhhi"|
                                   var == "shareactorscamp_all"|
                                   var == "rev_diffhhi.1"|
                                   var == "shareactorscamp_all.1") %>% 
  ggplot(aes(y = PE, x = var_num)) +
  geom_point(position = position_dodge(1/2) )+
  geom_errorbar(aes(ymin = lwr95, ymax = upr95),
                width = 0.05, size = .5, colour = "grey",
                position = position_dodge(1/2))+
  geom_errorbar(aes(ymin = lwr90, ymax = upr90),
                width = 0, size = 0.5,
                position = position_dodge(1/2)) +
  theme_bw() +
  geom_hline(yintercept = 0, lty = 2) +
  coord_flip() +
  scale_x_continuous(breaks = 3:4,
                     labels = c("Camp Diversity", "Camp Size")) +
  labs(x = NULL,y = "Standardized Estimate") +
  facet_wrap(~effect) +
  geom_text(data = data.frame(PE = c(0.16,0), var = spat_fx$var[c(3,4)],
                              var_num = c(4.25,4),ordering = c(2.09,2.09), 
                              lab = c(paste("rho==", round(mod1$rho,2),sep = ""), NA), 
                              effect = c("A: Direct Effects", 
                                         "B: Indirect Effects"))
            , aes(label = lab), size = 3, parse = T) +
  geom_text(data = data.frame(PE = c(0.16,0), var = spat_fx$var[c(3,4)],
                              var_num = c(4.15,4),ordering = c(2.09,2.09), 
                              lab = c(paste("(", round(rho_sd, digits = 2), ")",sep = ""), NA), 
                              effect = c("A: Direct Effects", 
                                         "B: Indirect Effects"))
            , aes(label = lab), size = 3, parse = T)
  


small_plot     

ggsave(plot = small_plot,
       filename = "CampResults.jpg", device = "jpg",
       width = 6.99, height = 4.99) 




